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Abstract —Probabilistic tractography based on diffu¬ 
sion weighted MRI has become a powerful approach for 
quantifying structural brain connectivities. In several 
works the similarity of probabilistic tractography and 
path integrals was already pointed out. This work 
investigates this connection more closely. For the so 
called Wiener process, a Gaussian random walker, the 
equivalence is worked out. We identify the source of the 
asymmetry of usual random walkers approaches and 
show that there is a proper symmetrization, which leads 
to a new symmetric connectivity measure. To compute 
this measure we will use the Fokker-Planck equation, 
which is an equivalent representation of a Wiener 
process in terms of a partial differential equation. In 
experiments we show that the proposed approach leads 
a symmetric and robust connectivity measure. 

1. Introduction 

Diffusion MRI has become a very important tool for 
understanding the living brain tissue ([E]). It can reveal 
both macro- and microstructural features of the neu¬ 
ronal network of the human brain. Tractography tries 
to characterize the structural connectome to understand 
the details of the interregional relationships of the hu¬ 
man brain. Tractography algorithms may be divided in 
deterministic, streamline-based methods m, m), proba¬ 
bilistic approaches m, [22] ) and global approaches ( [T9| , 
EH). In probabilistic tractography one basically draws 
samples from a distribution over paths and computes some 
statistics over these paths, for example, recording their 
endpoints. In some work the similarity to the notion of 
path integrals appearing in quantum mechanics |8| and 
statistical physics m was already pointed out m, 0, 
m, 0- Rigorous mathematical investigations show that 
the basic stochastic process behind path integrals is a 
so called Wiener process, a continuous Gaussian random 
walker [31]. In this work we will recap the foundations of 
the theory of Wiener processes and path integrals. Based 
on this we build a path integral that leads to a symmetric 
brain connectivity measure. We will see that the source 
of asymmetry of conventional walker principles is due to 
particle conservation. By dropping particle conservation 
and starting from the path integral perspective we will 
find a symmetric brain connectivity measure. Besides the 
walker perspective and path integrals, there is a third, 
equivalent approach, which describes expectation values 
of Wiener walkers by partial differential equations (PDF). 
This equivalence will be used to give an algorithm for 
the computation of the connectivity measure. It will be 


based on solving a large linear system, describing a mixed 
diffusion/convection process. We also introduce a novel 
discretization scheme to avoid the heavy directional de¬ 
pendency, which usually appears for discretized convection 
operators. PDFs in the context of diffusion MRI have 
been previously used for regularization [23], connectivity 
estimation [27] , [12] , [32] and fiber density estimation [25] . 
Our proposed PDF contains a convection operator, but in 
the joint position/orientation space, in [3] also convection 
was used for tractography, however just in position space. 
In [32] probabilistic tractography was put into a logical 
framework and and an algorithm was derived which is 
based on solving a large linear system. This approach also 
works just in position space such that crossings cannot 
handled adequately as long as only local neighborhoods 
are considered. The system matrix describes a anisotropic 
diffusion process. 

II. Method 

Apart from a few examples most methods to estimate 
brain connectivity are based on the walker principle. 
Fiber tracts are initiated from certain seed points and 
are iteratively built by following locally defined directions. 
While the deterministic tracking approaches are more used 
for illustrative purposes, the probabilistic ones are more 
related to quantitative connectivity analysis. The mathe¬ 
matical principle behind both approaches is an integration 
of streamlines along the underlying data. In probabilistic 
tracking the process can be seen as a Markov process. The 
iteration process is simple: if we call s(t) the current state 
of the tracker at step t, then the next state s(t+l) is drawn 
from some transition probability density IF(s(t + l)|s(t)), 
which may depend in some way on the DW-measurement. 
This process is basically a random walk in the state space 
of the tracker. If W is Gaussian, it is possible to formulate 
the limit for very small time steps, which results in a 
Wiener process, which we will now concentrate on. 

A. Wiener Proeess 

In physics and mathematics there exists a huge collec¬ 
tion of concepts for the analysis and characterization of 
Wiener random walkers. Basically, there are three perspec¬ 
tives which are mostly equivalent. 

The first one, which is closest to the above described 
Markov process is the concept of stochastic calculus and 
stochastic differential equations (SDF) [31]. In physics. 
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SDEs are usually written as Langevin equations, which do 
not have the strong mathematical footing of a SDE, but 
are simpler to read and more similar to ordinary calculus. 
To give an example (actually this example already covers 
everything we need later) we want to consider a simple 
diffusion process with some additional drift. Let the state 
of the ’tracker’ be s G and v a vector field causing the 
drift. Then, the corresponding Langevin equation is 

s{t) = v{s{t)) + T]{t) (1) 

where 77 stands for mean free white noise with variance 
per unit time. It is uncorrelated in time, meaning 
{rj{t)r]{t')) = — The dot in s means time differen¬ 

tiation. An approximative numerical integration gives rise 
to the following propagation scheme: 

s(t + At) = s{t) + v(s(t))At -|- u{t)^/At. ( 2 ) 

The first term v(s(t))At looks like ordinary Euler integra¬ 
tion. The second, stochastic term u(t) is drawn indepen¬ 
dently for each time step from Af{0^a‘^). The factor \/At 
expresses the fact that the process r]{t) is not differentiable: 
variances add up and not standard deviations. Note the 
difference between the continuous Wiener process r]{t) and 
the discrete process u(t) which is just defined at discrete 
time points kAt. Both are related by r]{t)dt = 

u(t)\/At. The transition probability W as described in 
the beginning for the discrete process described in (H is a 
Gaussian IT('|s(t)) = A/’(s(t) At v(s(t)), Atcr^). 

The second perspective describes expectation values of 
an ensemble S of random walkers described by 0 . Sup¬ 
pose we have generated such an ensemble S of walkers all 
starting at s( 0 ) = xq and we want to know the distribution 
of the states at some time T, i.e. 

Pt(x|xo) = y] <5(s(T) - x) (3) 

{s}e5 

Note, that the sum ranges over complete random walks or 
paths. To make this formally more transparent we always 
write curly brackets whenever we refer to a path (a chain 
of states of a walker) as a whole and not just a particular 
state at a specific time. In fact, the distribution described 
in 0 is a solution of a partial differential equation, which 
is called the Eokker-Planck equation (EPE) and is the 
master equation of the proposed continuous stochastic 
process. Eor this example, the EPE takes the form 

2 

Pt = -V- (vpt) + yApt = npt (4) 

where V is the usual gradient operator in and A the 
Laplacian. If we integrate this equation with respect to 
the initial condition po(x|xo) = S(x — xo) we just resemble 
the distribution given in 0 . The function pt{'x\xo) is 
also known as the propagator or Green’s function of 
the corresponding stochastic process. It can formally be 
written as pt = exp('Ht). Note that pt(x|xo) does not 
necessarily need to be normalized like a probability, e.g. 
walkers can die at the boundaries, which manifests in 
the boundary conditions of the corresponding EPE. The 


proposed algorithm will be based on discretized solutions 
of the steady state solutions of a symmetrized EPE. 

The third perspective is the path integral concept. Erom 
the viewpoint of a brain connectivity, the theory of path 
integrals is probably the most appealing one, however, 
they are probably the mathematical least understood con¬ 
cept and do not give a constructive way for designing an 
algorithm. However, it is essential for the understanding 
of what we are actually doing. The idea is to compute 
Pt(x|xo) by summing over all path starting at xq and end¬ 
ing at X weighted by its probability. Eor a rough derivation 
and motivation recall the discrete time representation. To 
compute the total probability that a particular path is 
drawn we take the product along the path 

N 

P(s(AAt),..., s(At)|s( 0 ) = xo) = IE(s(iAt)|s((i—l)At), 

but note that IT is a probability density. So, it is not 
totally correct to just multiply them, the problem will 
be discussed below. By taking the logarithm of the above 
product it can be converted into a sum and in the con¬ 
tinuous time limit N ^ 00 with NAt = T we have the 
integral form 

-log(P({s}|s(0) = Xo)) = - [ \og{W{s{t),s{t))dt 

Jo 

The functional L{s{t), s(t)) = — log(IT(s(t), s(t))) is some¬ 
times called the Onsager-Machlup functional (El, ID, 
na, p) or, in the style of mechanics, the Lagrangian. 

It describes the cost of a path. Let us further follow 
the naive way and use 0 and the corresponding W to 
compute L. By disregarding the normalization constant 
(which actually diverges in the limit), we can show that 

-^sym(s, s) = ~ ^(^)) (^) 

which seems to be plausible but is wrong. The correct 
answer is 

p(s,s) = ^(s-v(s))^ + i(V-v)(s) (6) 

The error is coming from the faulty assumption that we 
can simply multiply the conditional transition densities IT 
to get the total probability, which leads to the diverging 
normalization constant and the missing term V • v. Ac¬ 
tually we have to consider some small volumes instead 
of infinitesimal points. We have to consider the portion 
of walkers starting in a volume around s((i — l)At) and 
arriving in a volume around s{iAt). This leads to an 
additional factor (1 -|- which is the infinitesimal 

volume change during time At caused by the drift v. 

In more rigorous derivations of this result the factor is 
justified by the Jacobian of a change of variables (0, HI]), 
or by the considering the most probable tube around the 
path (P). 

Let us consider the final equation stating that the 
probability for a walker starting at s( 0 ) = xq and arriving 
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at s(T) = X is the sum over all paths connecting this two 
points weighted by their probabilities 


Pt(x|xo) = exp 


(-r 


L{s,s)dt (7) 


where denotes the set of all paths of length T 

starting in xq and terminating at x. In fact, the sum over 
all path is an integration in a functional space, for more 
details see e.g. i, m, ED- For convenience we omit to 
write here any normalization constant and assume that 
the measure is appropriately normalized. 


B. Reversibility and Symmetry 

In fact, the additional term (V • v)/2 in is quite 
important for us, because it is responsible for the fact 
that the described random walk is irreversible m)- In 
the context of brain connectivity this is directly related 
to the symmetry of the connectivity measure, so we will 
detail this out now. If we call V the Lagrangian with 
negative velocities v — v and s'(t) = s{T — t) the path 
traversed in opposite direction. Then, it is easy to see that 
J I/'(s', s')dt ^ f I/(s, s)dt with the L obtained in This 
also leads to pt(x|xo) 7 ^ p^(xo|x). This asymmetry is ba¬ 
sically caused by drift together with particle conservation. 
Considering the corresponding FPE in Q makes it more 
clear. While the diffusion term A is symmetric (selfadjoint, 
see Appendix IV-A for definition), the convection/drift is 
not antisymmetric, i.e. the operator B' corresponding to 
the flipped Lagrangian V is not equal to adjoint operator 
which has to be case to fulfill pt(x|xo) = p^(xo|x). 
However, for I/gym everything is different. In fact, we have 

/ -^sym(® ~ / Asym(s,s)(it 

Jo Jo 

which should also lead to a symmetric propagator 
Pt(x|xo) = p^(xo|x). But how to compute 0 for -Lgym 
in practice, what is the corresponding FPE? The term 
V • V acts like an additional potential energy. In fact, the 
Eeynman-Kac (EK) formula [17] tells us that such a poten¬ 
tial can be directly integrated into the EPE: suppose we 
have a given Lagrangian L together with its path integral 
and a corresponding EPE with operator B. Eor a potential 
field V{s) we define a new I/'(s,s) = I/(s,s) — E(s), 
then the EK formula tells us that the corresponding EPE 
operator is just B' = B BV . Applying this our problem 
with V = (V • v)/2 gives 


BsymP — Bp ■ 




v)p 


(vp) 


^(v-v)p-V 

y Ap - i(v ■ Vp + V ■ (vp)) 


Erom the last line we can also see the symmetry prop¬ 
erty of Bsym- The convection part V • (vp) + v • Vp is 
antisymmetric in the sense of the operator adjoint and 
hence the operator B^y^, which corresponds to flipped 


velocities, is just the adjoint of Bsym- One can also see 
that the particle conservation is lost: we cannot write 
Bsym as the a divergence of something. Note the similarity 
of Bsym with the quantum mechanical Hamiltonian of a 
particle in a vector potential, which is up to some constants 
(iV —v)^ = — A —^(Vv —vV)+v^. In quantum mechanics 
the Hamiltonian is indeed self-adjoint due to the additional 
complex unit in front of the anti-symmetric part. 


C. Steady States and Path Trails 

Up to now we have considered only paths with some 
specific length T. However, we are interested in all path 
without any length restriction. That is, we sum up over 
all paths connecting xq and x of arbitrary length to get a 
density p(x|xo) which is independent of T, i.e. 

nOO 

p(x|xo) = dT pt(x|xo) 

Jo 

= ^ g- !o L(s,s)dt 


where is the set of paths of any length connecting 

Xq with X. With the assumption hmT^ooPT(x|xo) = 0, 
i.e. all walkers eventually die, the function p(x|xo) is 
the steady state solution of the corresponding EPE, i.e. 
p(x|xo) is the solution of the equation 

-Hp(x|xo) = 6{x - Xo). 

This is the basic type of equation we will solve to esti¬ 
mate brain connectivities. In the following we call p(x|xo) 
the connectivity amplitude. Erom ordinary probabilistic 
tractography we know the so called length bias. With 
increasing euclidean distance, connectivity values decrease 
dramatically, usually exponentially. To account for this 
bias a linear reweighing was proposed in [4], that is, 

POO 

Piin(x|xo) = / dT T pt(x|xo). (8) 

Jo 

We will see later how pun can be computed in practice. 
Due to the exponential decrease one could also tend to 
use exponential corrections, i.e. use 

pOO 

P/t(x|xo) = dT e'^^pT(x|xo) (9) 

Jo 

with positive /i: G M. In fact, this results in a spectral shift 
of the operator B ^ B k. Of course, the choice of n is 
difficult. Large n will cause bad condition numbers for B 
and, in the extreme case, lead to a diverging p;^(x|xo). 

Until now, we were interested in a measure how strong 
two points xi and xq are connected, but one might also be 
interested in how the actual connection looks like. Instead 
of just summing over all paths connecting two points, one 
can additionally image the path itself by counting how 
often the path visited a specific position x. We call this 
image the trail image of a path s, and it is defined by 

rs(x) = [ 6{-k — s{t))dt 
Jo 
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Now, we can write the mean trail image of a path connect¬ 
ing two point xi and X 2 by the expectation value 

In order to compute f we have to use the so called Einstein- 
Smoluchowski-Kolmogorov-Chapman (ESKC) relation, or 
simply semigroup property, which tells that the path- 
integral, or the corresponding propagator can always be 
split like 

Pt(xi|xo) = j dx.' PT-t'(xi|x')Pt'(x'|xo) 

for any given intermediate point 0 < t' < T. This means 
we can slice the propagator at any given length t' and 
consider all possible intermediate positions s(t') = x' the 
path has at f and integrate over them by accounting for 
the probability of the path to get from xq to x' and to 
reach the target xi starting from x'. To compute f we 
first consider the point density of the path at some specific 
length t'. That is, we count how often a path traverses the 
point X at length t' by 

pT,t'{^) = y^^x'(5(x-x')pT-t'(xi|x')pt'(x'|xo) 

= pr-t'(xi|x)pt'(x|xo) 

If we integrate now over all intermediate times 0 < 
t' <T and, again, integrate this over all possible 0 < T < 
oc we get the mean trail image 

OO nT 

dT / dt'pT,t' (x) 

Jo 

Plugging this altogether gives the simple final result 


joint space of position and orientation s = (r,n), where 
r G and n G 5'2- The data we get from the diffusion 
measurement is basically an fiber orientation distribution 
/(r, n) defined on this joint position/orientation space. In 
fact, in this joint space there are several ways to define 
meaningful connectivity measures. We propose to use the 
following Lagrangian 

L(r, r, n, h) = ^(r - n/(r, n))^ + (10) 

as path-costs. Paths with s ^ n/ and small changes in 
n are assigned with small costs, and hence, with high 
probability. Erom the walker perspective, there is a con¬ 
vective force which drives a walker with internal state 
n into direction n, while the speed is proportional to 
the data term /. Additionally there is diffusion on the 
orientation variable, enabling the walker to change its 
directions, or conversely, penalizes too strong bending. 
Imagine the system as a network of pipes oriented in all 
possible directions carrying a fluid which flows with speed 
/(r, n). The flow is not perfectly convective, neighboring 
parallel pipes can exchange fluid (this is spatial diffusion) 
and also crossing pipes exchange some fluid (this is the 
angular diffusion). In fact, we can also let spatial diffusion 
ar ^ 0 such that there is only pure convection. However, 
as we will see later, the discretization of the EPE will force 
us to have finite values for cr^. 

The EPE operator is not totally straight-forward, be¬ 
cause there is an additional term coming from the curva¬ 
ture of the sphere m, ID)- It looks like 

An) - i(v ■ Vr + Vr ■ v) + ^ 


Txi|x2 


)(x) = y 

Jo 


T(xi|xo)(x) =p(xi|x)p(x|xo) 

and, as we already know how to compute p we are ready. 
Eurther note that 

pOO 

c^x f(xi|xo)(x) = / dTTp(xilxo) 

Jo 

= Plin(xilxo), 

that is, if we integrate the mean trail image, we just get 
the linearly reweighed version of (§ and have also a way 
to compute ^ in practice. 

D. Application to DW-MRI data 

All concepts were introduced by considering a 3D dif¬ 
fusion process with drift. But, how to apply this kind 
of stochastic process for brain connectivity analysis? The 
data is basically tensor-like, meaning it is point symmetric, 
there is no velocity field available. Even when we have 
orientation distributions, there is still no convection-like 
force. So, what to do? 

Eirst of all, we have to realize that usual probabilistic 
tracking algorithms are not Markovian with respect to 
the position r G The propagation direction usually 
depends not just on the position, but also on the previous 
step made. So, the state space of the tracker should be the 


where the velocity is defined by v(r,n) = n/(r, n). The 
additional constant ^ is coming from the curvedness of 
the sphere which is i? = 2/p^ for a 2-sphere of radius p. 
However, there is actually no metric connection between 
the position space of the orientation space 82 - The 
’radius’ of the sphere 82 has no actual meaning. So, the 
value of R/12 is rather arbitrary. In fact, R/12 acts like 
an additional exponential weighting, like the k we already 
introduced above and could be absorbed by that. So, we 
exclude R/12 in the following, or, imagine p to be pretty 
large. 

Let us now look at the symmetry properties of %. Let Z 
the operator that does a point reflection on some function 
a(r, n) in the sense 

(^a)(r,n) := a(r,-n) 

Then, we can easily show that 

nZ = 

because the data / is symmetric Zf = /. In other words, 
the operator RZ is selfadjoint. Eormally we can write the 
connectivity amplitudes as 

p(r,n|ro,no) = -(e(r,n) |^“^e(ro,no)) 
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where e(r/,n') (^5 — T)6{n' — n). For the connec¬ 

tivity amplitude the symmetry means 

p(r, n|ro, no) = p{ro, -no|r, -n), 

because = Z{'H'^)~^Z. 

If we now think of some function a(r, n) describing a 
seed point density, and 6(r, n) a terminal point density. 
For both we can think of indicator functions describing 
some cortical region of interests. Then, we can define the 
symmetric connectivity measure c(a, h) to be 

c(a,6) = -{a\Zn-^\h) 

which can be written in terms of the connectivity ampli¬ 
tudes p(r, n|ro,no) as 

c(a, b) = j j a(r, -n)p(r, n|ro, no)6(ro, no) 

So, c{a,b) is the path integral over all path connecting 
region a with region b according to the Lagrangian L. 
Accordingly, we can also define c'^ and corresponding 
to piin in equation (§ and in equation respectively. 
The path trail image is 

T'(a|6)(i’,n) = (a|'H“^er,„)(er,„|'H“^6) 

= {er,n\2n-^Za){er,n\n-^b) 

= {er-n\'H~^Za){er,n\'H~^b) 

Note that '^{Zb\Za)' 

Usually the seed functions a and b do not depend on 
n because we do not have any preference for the starting 
orientation, so Za = a and Zb = b. Several relationships 
are then simplified, e.g. Zf(^a\b) = U6|a) c{a,b) = 
{a\'H~^\b). For the special case if aro(i*,n) = (5(r — tq) 
we can also define the connectivity measure between two 
spatial locations r and fq as 


c(ro,ri) = c(aro,ari) (H) 

In all experiments we will report the normalized connec¬ 
tivity amplitudes 


Cn{a,b) = 


c(a, b) 


y/c{a,a)c{b,b) 


which is quite natural and reminds one of Pearson’s cor¬ 
relation coefficient. 


E. Angular Constraints, the Data Term and Boundary 
Conditions 


The Lagrangian (10) penalizes bending, which is a good 
prior for straight fibers, but underestimates curved fibers. 
In [4] this problem does not appear, because the walker is 
always back projected onto the nearest fiber direction in 
a voxel. Such a behavior cannot be realized by a Wiener 
process. However, we can mimic such a behavior by an 
angular constraint such that the walkers do not deviate 
too strong from the main fiber directions. That is, we only 
simulate paths where n is not too far away from the main 
fiber directions, which is similar to the maximum angle 
thresholds known from ordinary streamline algorithms. 


but not so rigorous (it is still possible to find paths where 
the spatial tangent r to the path is not perfectly n). 
Now, we are able to set pretty large, i.e. low penalty 
on bending, but the walker will be along the main fiber 
directions. In fact, this idea will also help to keep the 
running time and memory consumption of our algorithm 
in a reasonable range, because we only have to simulate the 
domain in the neighborhood of the fiber directions and not 
on the complete sphere. Note that the angular constraint 
may also be seen as an additional cost in the Lagrangian, 
i.e. the forbidden path are assigned with infinite costs. 

There are lots of ways to determine the speed function 
/, which represents the DW-data. For example, we could 
directly use the fiber orientation distributions (FOD), e.g. 
estimated by spherical deconvolution [28]. Or, one could 
estimate the main fiber directions by a local maxima 
detection of an FOD, and use them as anchor directions to 
construct the speed function. We opt for latter, primarily 
because the angular constraints can be much beter con¬ 
trolled. Let di be the local maximas of the FOD. Then, 
we construct the speed function by 

N 

/(“) = 

i=l 

with some fixed n > 0. The simulation domain is now 
defined by thresholding this speed function, i.e. we only 
consider regions with / > e for some e > 0. 

To solve the FPE we have to set some boundary condi¬ 
tions. It is natural to let the walkers die once they touch 
the boundary of the domain, which is equivalent to a 
zero Dirichlet conditions at the boundary. So the complete 
problem we want to solve is 

— Bp = a where p{dfl) = 0 (13) 

where dft denotes the boundary of the simulation domain, 
p the unknown steady state solution and a the seed region 
where the walkers are emitted. Of course, other boundary 
condition, like Neumann conditions, might be applied, 
for example at the transitions to the ventricles. However, 
we restrict here all considerations to the simple Dirichlet 
conditions. 

F. Implementation, Diseretization and Solver 

In order to solve the equation —TLa = b we have to 
discritize the operator TL. The convective part of TL makes 
this a pretty hard problem, usual Finite Element Methods 
are only applicable for diffusion dominated problems. So, 
we decided to use a finite difference upwind-downwind 
scheme. However, such schemes are known to introduce 
errors (known as false diffusion) when the convection 
direction is not aligned with the underlying discretization 
grid. In particular for steady state solutions these problems 
become severe and the solutions show heavy orientation 
dependency. Eortunately our problem is special, the con¬ 
vection direction is always fixed for a volume p(r, n) 
with fixed n, namely n. So, we propose to steer the grid 
for each direction along the current convection direction. 
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Practically, we discretize the sphere into N direction 
homogeneously distributed over the full sphere (in the 
experiments we used N = 128 directions). So, we have 
N volumes, where the discretization grid of the volume 
associated with the ith direction is steered such that 
one of the coordinate axis is along direction . In this way 
we are able to produce quite clean convection, however the 
angular diffusion gets disturbed. A voxel at some discrete 
coordinate in volume i does not have a unique partner in 
volume j at exactly the same discrete position. We need 
this partner in order to implement the spherical Laplace- 
Beltrami operator An, which only acts on the spherical 
coordinate. The simplest way to account for that is to 
use trilinear interpolation, which, unfortunately causes 
another unwanted spatial diffusion effect. However, this 
false diffusion is not so disturbing: it is independent of the 
orientation and the error is proportional to the angular 
diffusion and not to the speed function. Second, we can 
reduce the error by a finer discretization of the spatial grid. 
We quantified this additional spatial diffusion in an exper¬ 
iment and found that a angular displacement of tt = 180° 
(for a A = 128 points on the sphere) causes an additional 
spatial displacement of about 5 spatial grid units, i.e. the 
cost of a 180° turn on the spot compares to a 5 voxel jump. 
So, for some given theoretical and cr^, the actual 

turns out to be approximately = ^/cr^ + (dnS/Tr)^, if 
ar is given in grid units. As this is already pretty much we 
decided to set the theoretical to zero in all experiments. 
In general, for arbitrary number of p oints on the sphere , we 

found the approximation — \/ (<^nV^/7)^, i.e. 
with a finer sphere discretization the effects get stronger. 
But again, note that the ’physical’ amount of this false 
diffusion depends also on the spatial grid size. All the 
gritty details about the discretization scheme are given in 
Appendix |IV-B 

The whole matrix is scattered with MATLAB’s sparse 


matrix capabilities. To solve the discretized equation (13) 


we found the GMRES algorithm the most stable and effi¬ 
cient |26], we just used the MATLAB’s gmres command. 
As number of restarts we used 5. The tolerance value is set 
to 10“^ in all experiments. The running time depends on 
the size of the problem. If we use the setting as described 
below in the experiments we get usually a system with 
about 2—3-10^ variables, which is scattered in 1—2 minutes 
on a common Desktop PC (Intel 17, 16GB). Solving the 
equation also takes about one minute. 


III. Experiments 

We want to begin with a small demonstration of the 
proposed discretization scheme. Let d some arbitrary di¬ 
rection and the speed function is constructed according to 
([T^ with n = 25 and e = 0.02. In Eigure[^we choose d = 
{cos{ip), sin{ip)^ 0) with ip = 0.37r (a,b,e) and p = tt (c,d,f). 
Eurther we choose = 0 and an = 7r/12. Eigure[^ (a-d) 


shows the solution to (13) with a(r, n) = (5(r — tq) where fq 
is located in the lower left corner of the simulation domain. 
The spatial intensity maps in Eigure[^show the solution 



Fig. 1. Demonstration of the proposed discretization scheme. In (a- 
d) the novel scheme for different granularities are shown. In (e-f) the 
ordinary scheme is shown for comparison. 


p integrated over the angular variable f^^p(r,n)dn, i.e. 
we actually show the connectivity amplitudes c(r, fq) as 
defined in 0- The matrix size is of size 30 x 30 x 7. In 
Eigure[^,b,c) we used a sphere with N=128 directions in 
(d) with N=256 directions. Eurther, we demonstrate the 
effect of spatial oversampling, i.e. a upsampling by a factor 
of 2 means that the actual matrix size is in the range of 
60 X 60 X 14 (for details see Appendix IV-B). By comparing 
(a) and (b) we can obviously see how the the upsampling 
reduces false diffusion. Similarly comparing (b), (c) and 
(d), we see that the amount of false diffusion does not 
depend on the direction, but the directions are not perfect 
due to the discretization of the sphere. Eor comparison we 
also simulated an ordinary up/downwind scheme in (e),(f) 
with an = a^ = 0. The amount of false diffusion depends 
dramatically on the direction. 


A. Numerical Phantom 

Eor further demonstration we constructed a numerical 
phantom consisting of crossing and bending configurations 
(see Eigurej^)). Six seed locations fi, ..., Fe were selected 
such that they are pairwise connected as (1-2), (3-4) and 
(5-6). We did not simulate the MR-signal, but created 
the underlying directions directly. The directions are 
created continuously (not aligned with the discretized 
sphere directions), and a pseudo EOD is generated using 
exp (a ((dill) ^ — 1), which is shown in Eigurej^). Eor 
robustness analysis we distorted the directions by Gaus¬ 
sian noise of standard deviation a^z- 

Eigure c) gives some first results of our approach. 
We show the connectivity amplitudes c(F,Fi) as images, 
where f^ is one of the fixed seed point. Eigure |^) shows 
c without noise and b) with noise {a^z = 0.1). Eigure [^) 
shows the same experiment but with exponential length 
bias correction of k = 0.01. Examples for the path trails 
images are shown in Eigure [^) that correspond to the 
experiment in Eigure [^). Trail images for all pairs of 
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Fig. 3. Quantitative analysis of connectivity values (c-values) over 50 trials. The barplots show mean c-values together with standard 
deviation. The x-axis indicate all possible seed pairs (i,j) for z, j = 1,..., 6 of the numerical phantom. 



d) Phantom, seed locations e) Path trails 


Fig. 2. In (a-c) connectivity amplitudes for three different settings 
are shown, the color scaling is fixed and the same for all three 
settings. In (d) the FOD of the phantom together with the seed 
locations is given, and (e) shows the matrix of all possible path trails 
together with the connectivity amplitudes in logarithmic scaling in 
the diagonal. 


seed combinations are shown, the diagonal shows again 
the connectivity amplitudes but in logarithmic scaling. 


To get a quantitative picture we repeated the above ex¬ 
periment for different noise levels = 0.1, 0.2 and differ¬ 
ent discretization scheme, namely spherical discretization 
with N = 128 and N = 256 directions and with and w/o 
spatial upsampling (with a factor of 2). Additionally, the 
orientation of the discrete sphere was randomly changed 
for each run. As reference approach we followed [4] (for 
details see Appendix IV-C), with As = l,cr = 0.2, 
without revisits and with length bias correction. Figure 


shows barplots of the normalized connectivity amplitudes 
Cn(ri,rj) for all pairs of seeds averaged over 50 runs 
together with the standard deviation. 


B. In- Vivo Human Brain 

To investigate our approach on real DTI data we 
considered 28 scans of healthy volunteers at a b- 
value of 1 ms/linn? with 61 diffusion directions and an 
isotropic resolution of 2 mn?. A white matter prob¬ 
ability map was generated with SPM (Version 8, 
http://www.fiLion.ucLac.uk/spm/, [10]) on a Tl-weighted 
scan, which was co-registered to the 6o-scan of the diffusion 
sequence. For each subject the scans were repeated two 
times (in two different sessions) to allow to investigate 
the robustness of the approach. To compute the FOD 
a constrained spherical deconvolution was used together 
with a I/i-regularizer m- The fiber response function 
(FRF) was chosen as FRF{n) = ex.p{—bDi{nf^ni/) with 
Di = liiw?/ms for all subjects and reconstructed on a 
V = 128 sphere. The local maximas of the FOD are 
determined by fitting a quadratic form to the neighbor¬ 
hood of the discrete local maxima (the same approach 
as used in ESI). Spurious local maximas were filtered out 
by neglecting all maximas smaller than 20% of the global 
one. The speed function is constructed in the same way as 
for the numerical phantom. The spherical discretization is 
kept at V = 128 and no spatial up sampling is applied. 
The probabilistic white matter segmentation of SPM was 
thresholded at 0.5 and used as white matter mask. To 
compare and analyze the group, the AAL atlas m was 
registered to the subjects’ native space, normalization was 
done with SPM 8. The first 90 ROIs (distributed over 
the neocortex) of the AAL atlas were used to compute 
connectivity matrices. 

Figure |^-e) show examples of the mean connectivity 
matrices (CM) over the whole group. The first 45 ROIs 
belong to the right hemisphere, the last 45 the right 
hemisphere. The CM is shown in logarithmic scaling. To 
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Fig. 4. In (a-e) average connectivty matrices (CM) for different approaches are shown. The matrices in (f,g) show ICC values for individual 
c-values. Plot (h,i) display the ICC for all c-values (all ROI pairs) in the CM. Plot (j) shows agreement values. 


get a comparable contrast we used the following formula 
ciog(a,6) = log(t + c^(a,6)) where t is the 20% quantile 
of Cn over all regions. The CMs obtained by our approach 
are compared to the reference approach [1] (Figure|4ji-e)). 
To investigate the robustness we computed the intraclass 
correlation coefficient (ICC]0 for each connectivity value 
(c-value) in the CM. In Figure [^) we show the ICC 
over all c-values in the CM sorted by ICC magnitude 
for different settings: we varied dn and show one result 
with exponential weighting = 0.01). For all results we 
also show the corresponding linear weighted result 
Figure [^) compares our approach with [4] with different 
parameters: a = 0.2, 0.5, As = 1, with (w/o) revisits and 
length bias correction (Ibc). In Figure|^,g) we show for our 
approach and the reference the ICCs for the individual c- 
values as a matrix. Table U shows mean and median ICCs 
for our approach and the reference with more different 
settings. It is common to put thresholds to get a binary 
decision of connectivity. To investigate the robustness 
of such a thresholding operation an agreement measure 
is calculated as follows: let t be some threshold, then 
a = I Cl <t/\C 2 <t\ denotes the number of regions that are 
not connected in scan 1 and in scan 2, where the number 
is counted over all possible ROI pairs and subjects. The 
agreement value a is normalized by the total number of 
non-connected regions s = \ci <t| + |c 2 <t|, and hence, 
a/s is a number between 0 and 1. If we now vary over all 
thresholds t and plot a/s as a function of s we get plots 
displayed in Figure]^). 

Finally, Figure shows examples of path trail images for 
two ROI pairs of the AAL atlas. The two ROIs are high¬ 
lighted in green and blue, respectively. The white matter 
mask is depicted in dark gray. The path trail r(a| 5 )(r, n) 

^If Cl and C 2 are c-values in the CM for scan I and scan 2, respec¬ 
tively, then the ICC is ICC{c) = {(ci—c)(c 2 —c)}/{(c—c)^}), where {} 
denotes the expectation over the whole group and c = {(ci + C 2 )/ 2 }. 
The ICC is I if all the variance of the c-value can be explained by 
differences amongst the subjects. 


TABLE I 

Quantitative Comparison by ICC 


Method 

Mean 

Median 

Our Approach 


an = 0.2 

0.47 

0.51 

n = 0.5 

0.64 

0.69 

an = 1.25 

0.73 

0.78 

an = 0.5, K, = O.OI 

0.67 

0.72 

an = 0.2, lin. weighted 

0.50 

0.54 

an = 0.5, lin. weighted 

0.67 

0.71 

an = 1.25, lin. weighted 

0.73 

0.78 

an = 0.5, R = 0.01, lin. weighted 

0.69 

0.73 

Reference [4] 


a = 0.2, As = 1 

0.65 

0.68 

a = 0.5, As = I 

0.55 

0.63 

a = 0.2, As = I, w/o revisits 

0.61 

0.64 

a = 0.2, As = I, no Ibc 

0.64 

0.67 

a = 0.2, As = 0.25 

0.55 

0.58 

0 - = 0.2, As = 0.5 

0.60 

0.63 



Fig. 5. Examples of mean path trails for two ROI pairs in the frontal 
regions. Green and Blue depict the ROIs. The path trails are shown 
as glyphs and averaged over directions. 


is obviously a function in the joint position/orientation 
space, which is shown in glyph representation, but note 
that the glyphs are not pointsymmtric as one used to 
observe. The traversal direction is in this example from 
the blue to the green ROI. 

IV. Discussion 

We proposed to use a Wiener process for the estima¬ 
tion of structural brain connectivities. The path integral 
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perspective gave us the insight why the usual random 
walker principle leads to asymmetric connectivities. To 
symmetrize the solution an additional potential was in¬ 
troduced, which symmetrized the path integral and the 
corresponding Fokker Planck equation. However, particle 
conservation and, hence, the usual walker perspective is 
lost. To solve the FPE we introduced a novel discretization 
scheme which avoids the directional dependency of the 
false diffusion introduced by the discretization of the 
convection operator. 

Initial experiments (in Figure show that the new 
discretization approach is valid. In Figure the visual 
results on the numerical phantom show the approach 
works also nicely under noisy conditions and that the expo¬ 
nential length bias correction is reasonable. Quantitative 
experiments show that it leads to a robust connectivity 
measure. Compared to [4] the new approach shows smaller 
connectivity values for bending configurations, which is 
due to the bending penalty which is inherent to the 
new approach (see Figure]^. One can also see that the 
robustness of our approach does not really depend too 
much on the discretization granularity. Compared to [4] 
it is in most cases higher (lower standard deviation over 
the trials). For stronger noise {a^z = 0-2) our approach 
gives relatively higher c-values to false connections, but 
the c-values itself are more stable. Another difference to 
[4] is the way the curvature is controlled. If we want 
for our approach a higher connectivity for bending fibers 
(ROIs (5,6)), then we also get higher c-values for (1,3) and 
(2,4). For [4] it is different, the step directions are always 
back projected onto the nearest most likely fiber direction, 
which results in higher c-values for the bending (5,6). 

To investigate the new approach in-vivo, we considered 
a group of volunteers which were all scanned twice. We did 
a ROI-based approach based on the AAL-atlas. Overall, 
the CMs obtained from our approach and the reference 
are similar. As one might expect, with our approach strong 
bending fibers (like the connection between ’Occipital Sup 
L’ to ’Occipital Sup R’) are assigned smaller connectivity 
values. Analysis of the intraclass correlations of c-values 
shows that the new approach is quite robust. Usually 
higher angular diffusion gives more robust c-values. 
Robustness is further increased by a length bias correction, 
either linear or exponentially. Compared to [4] our ap¬ 
proach shows good performance. For [4] increasing angular 
diffusion does not help so much. One can see that in 
regions with low connectivity [4] has problems (see Figure 
|^,g)), which goes inline with our finding in the numerical 
phantom. We found [4] to work best for As = 1, a = 0.2, 
revisits are allowed and with length bias correction (see 
Table 1^. In this setting the median ICC is 0.68 which is 
acceptable, 50% of the c-values have a ICC higher than 
0.68. The highest value we can achieve with our approach 
is a median of 0.78. The agreement scores obtained from 
thresholded c-values show a similar but more pronounced 
picture (see Figure]^)). From Figure]^ it can also be seen 
that the estimated path trail images are meaningful and 
follow the expected pattern. 


Appendix 

A. Operator Formalism 

The space of square integrable function on is a 
Hilbert space H with the inner-product 

{q\p) = f 9 *(x)p(x)dx 

Sometimes we adopt Dirac’s bra-ket notation, i.e. by \p) we 
denote a vector in i7, by {p\ a linear 1-form. The adjoint 
of a linear operator A is defined to be the operator 
such that {q\Ap) = {A'^q\p) for all functions q^p e H. 


B. Discretization 


Let m(r) some predefined white matter mask. For each 
direction on the sphere, we create a voxel grid rU 
such that the x-axis is aligned with direction n^. Only valid 
triples (x, z) G are considered where ^ is on. 

The voxel spacing h = y z) ~ ^Ix y z) I chosen 

arbitrarily. In the experiments we used h = 2mm which 
is resolution of the DTI scan, and h = 1mm, that is an 
oversampling factor of 2. The sphere has a natural neigh¬ 
borhood system Af{i) obtained by Voronoi tessellation. 
Let Pi^x,y,z be some function on the introduced discrete 
domain. To discretize the spherical Laplace-Beltrami 
operator we define the discrete operator A^ as 


{^nP)i 


E 

keM{i) 


Pi,: 


- E 


'^k,a,b,c Pk,a,b,c 


(a,b,c)e 


where is neighborhood system needed for tri- 

linear interpolation at position ^ in the k-th volume, 

and Wk^a^b,c are the corresponding weights. However, the 

matrix A^ is not symmetric in the sense of the matrix 

transpose due to the one-sided interpolation. So the final 

2 ~ ~ 

discretized spherical Laplacian is ^(A^ + A^). The factor 
A is empirically determined to get standard deviations of 
approximately <7^. 

To implement the convection part, the speed function 
/(r, n) given in the native DTI frame is also transferred 
by trilinear interpolation onto the new coordinate frame 
fi,x,y,z = f{'^lxyz)^'^'i')' "^he simulation domain is deter¬ 
mined by thresholding fi^x,y,z' For the convection operator 
we have to distinguish between ’real’ boundaries and 
boundaries coming from the thresholding. The discretized 
convection generator V is defined as 

{ Pi,.-i,y,z-Pi,.,y,z otherwise 

Pi,x-l,y,z/h if hi^x+l,y,z = 1 

-P^,x,y,z/h Ah^x-l,y,z = -l 

Pi,x,y,zIh if bi^x — l,y,z — 1 

where bi^x,y,z is an boundary indicator, where 6=1 means 
the boundary is coming from thresholding and 6 = — 1 
from a ’real’ boundary. Finally, the convection part of the 
equation is — (FV -1- VF), where F is the matrix with 
fi,x,y,z on the diagonal. 
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C. Reference Method 

For comparison we followed [4] . Suppose for some given 
fiber directions, which are always extracted in the same 
way as for our approach. The state of a walker is its current 
position r and its last step direction n. The propagation 
law for a single walker is: The current voxel contains a 
set of fiber directions d^. Determine the fiber direction d 
which is closest to the current n. If the angle between those 
is above a threshold of 80° tracking is stopped. Otherwise, 
disturb this fiber direction by Gaussian noise d^ = d + 
rjcr with deviation a. Renormalize d^ to unity and track 
along this direction with step width As, i.e. = r + 
Asdn and set = d^. If the new position is out of 

the tracking mask, stop tracking. If desired, one can also 
stop tracking whenever a voxel is visited more than once. 
During tracking we count how often a walker has visited a 
voxel resulting in a probabilistic map (PM). Alternatively, 
to account for the length bias, one can also weight with 
visitation count by the mean length of all walkers visited 
the voxel. To compute the connectivity value between ROI 
A and ROI B, N walkers are ejected in each voxel in ROI 
A. The connectivity value is then the sum over the PM in 
ROI B. To get a symmetric value the connectivity from A 
to B is averaged with the connectivity from B to A. To set 
N we computed the intra-scan ICC, i.e. for each subject 
and scan in the group we started the algorithm twice and 
computed the ICC between the two algorithm runs. We set 
A, such that the median ICC over all ROI pairs is above 
0.95, which resulted in N = 5000 per voxel. 
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